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Abstract 

Power systems are undergoing unprecedented transformations with the incorporation of larger amounts of renewable energy 
sources, distributed generation and demand response. All these changes, while potentially making power grids more responsive, 
efficient and resilient, also pose significant implementation challenges. In particular, operating the new power grid will require 
new tools and algorithms capable of predicting if the current state of the system is operationally safe. In this paper we study 
and generalize the so-called energy function as a tool to design algorithms to test if a high-voltage power transmission system is 
within the allowed operational limits. In the past the energy function technique was utilized primarily to access the power system 
transient stability. In this manuscript, we take a new look at energy functions and focus on an aspect that has previously received 
little attention: Convexity. We characterize the domain of voltage magnitudes and phases within which the energy function is 
convex. We show that the domain of the energy function convexity is sufficiently large to include most operationally relevant and 
practically interesting cases. We show how the energy function convexity can be used to analyze power flow equations, e.g. to 
certify solution uniqueness or non-existence within the domain of convexity. This and other useful features of the generalized 
energy function are described and illustrated on IEEE 14 and 118 bus models. 


I. Introduction 

Power systems are experiencing revolutionary changes. Integration of renewable generation, distributed generation, smart 
metering, direct or price-based load-control capabilities have contributed to the revolution, but are also posing significant 
operational challenges by making the power system inherently stochastic and inhomogeneous. With these changes, the system 
operators will no longer have the luxury of large positive and negative reserves. Moreover, operating the future power grid 
will require developing new computational tools that can assess the system state and its operational margins more accurately 
and efficiently than current approaches. Specifically, these new techniques need to go beyond linearized methods of analysis 
and ensure that the power system is safe even in the presence of large disturbances and uncertainty. In this paper, we study the 
energy function approach and argue that this classic tool can be generalized to help solve a number of existing and emerging 
problems in operating modern power systems. 
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Introduced more than 50 years ago by Aylett | Ayl58| , energy functions have been used to develop direct methods of the 
power systems stability and analysis. Essentially, energy functions are physically motivated Lyapunov functions for certain 
simplified models of the power system dynamics. In the 70s and 80s, many results were derived for stability analysis of power 
systems,first via approximate Kron-reduction techniques and then using more accurate structure-preserving energy functions 
IIBH81I . Developments of this period are summarized in 0VWC85I . IIPai89l . Since then, much of the work has focussed on 
developing effective heuristics to produce non-conservative estimates of the region of power system stability (see IlChilll and 
references therein). Interest in these approaches declined in the 1990s and 2000s, partially due to the difficulty of generalizing 
these approaches to more detailed power system models, but also due to advances in computing technologies that made stability 
analysis based on real-time dynamic simulations feasible. 

Our results are also linked to more recent lines of research related to (a) necessary and sufficient conditions guarantee¬ 
ing existence and uniqueness power flow solutions IIIli92IIIASV81lllLSP99IIIDB14IIIDB12ll||DCB13l ; (b) analysis of the fixed 
point characterizations of the power flow equations and analytical approximation schemes for power flow solutions MBZ141 : 
(c) semidefinite relaxations of optimal power flow IILowl4bl . IILowl4al| and various extensions (see for example BLS B1 41 . 
I1MAL14I ). We argue in this manuscript that energy functions can be viewed as a unifying tool for several operational problems 
in power systems. This is based on the fact that for lossless systems (resistance of transmission lines are ignored) stationary 
points of energy functions are solutions of the Power Flow (PF) equations. Hence, many traditional operational problems 
dependent on the notion of the (static) PF equations, such as the PF analysis itself, the Critical Loadability Problem (CLP) and 
Optimum Power Flow (OPF) problem,can be studied using the energy function methodology. This viewpoint, first presented in 
EVU, I1BBC13I . is extended here both in terms of theory and applications. Moreover, we show that this viewpoint extends to 
the more general cases where the energy-like functions are not true Lyapunov functions, but simply functions whose stationary 
points correspond to PF solutions. This extension will allow us to deal with lossy systems that have fixed resistance to reactance 
ratio for all transmission lines. 

In order to obtain efficient formulations, we study a property of the energy functions that, to the best of our knowledge, 
has previously received little attention: convexity. Since power systems are designed to operate at a stable equilibrium point, 
the energy (Lyapunov) functions characterizing their behavior are obviously convex at least in a small vicinity of the stable 
equilibrium. However, to the best of our knowledge, the exact domain of convexity has never been characterized. The classical 
papers by Bergen and Hill IIBH81I on the structure preserving model of power systems, extended in flNM841 . IICRP85II to 
account for effects of the reactive power flows and variable voltages, form a solid base for our analysis. We study various 
assumptions under which the energy function is provably convex. It is well-known that when all nodes are the so-called 
(P, V) nodes (active power consumption/production and voltages are fixed/known), the energy function is convex 0BBC13I . 
flBV12ll as long as all phase differences are smaller than 90 degrees. In this paper, we generalize this result to the networks 
with {P, Q) nodes (i.e. nodes where active and reactive power injections/consumptions are fixed while voltage magnitudes 


are variable). 

Our main technical result is a description of the domain over which the energy function is convex. We provide a nonlinear 
but convex matrix inequality condition on the voltage phases and magnitudes that guarantees convexity of the energy function. 
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Moreover, we show that this condition provides an exact characterization of the domain of convexity of the energy function for 
the case of acyclic networks. For mesh networks, in general, the matrix inequality condition provides an inner approximation 
of the convexity domain. The condition for convexity is quite natural and captures the intuition that if voltage magnitudes and 
voltage phases at neighboring nodes are “not too far” from each other then the underlying energy function is convex. 

The rest of this paper is organized as follows. In section [I^ we describe the mathematical background and power system 
models. In section m we first analyze convexity of the energy function jointly in voltage magnitudes and phases, and then 
describe how the convexity can be used to design a provably convergent algorithm solving the PF equations (Section [lII-C| l. In 
section ??, we discuss related work and describe the advantages of the energy function methodology. In section |V] we provide 
numerical illustrations of our approach on the IEEE 14 and 118 bus networks. Section |Vl] contains conclusions and a brief 
discussion of the future work. Finally, Appendices |VIII-A| and |VIII-B| are reserved for the proof of the main theorem and a 
generalization to a special class of lossy networks. 

II. Modeling Power Systems 

In order to guarantee the existence of an energy function, we will need to make the standard assumptions used in the 
structure preserving models IIBH81I . which are fairly well justified for the high-voltage transmission network; 

(1) All transmission lines are purely reactive (network is lossless), i.e. the network resistive power losses are considered small 
and thus ignored. 

(2) Nodal active and reactive power injections/consumptions are constant. 


A. Notation 


The transmission network is modeled as a graph (V,£) where V is the set of nodes and S is the set of edges. In power 
systems parlance, the nodes are called buses and the edges are called lines (transmission lines). We shall use these terms 
interchangeably in this manuscript. Nodes are denoted by indices i = l,...,n = |V| and edges by ordered pairs of nodes 
We pick an arbitrary orientation for each edge, so that for an edge between i and j, only one of {i,j) and {j,i) is in 
S. If there is an edge between buses i and j, we write i ^ j,j ^ i. Sometimes, it will be convenient to denote edges by a 
single index k. In this notation, ki and ^2 denote the “from” and “to” ends of the edge k. Define the edge-incidence matrix 
A G as follows: 


1 if A: S fi G V, fci = i 


^ik — 


-1 if fc G £J G V,fc2 = j G V 


0 otherwise 


The set of vertices includes all load buses as well as buses representing generators (a single bus, or terminal and internal 
bus per generator). This allows us to develop a unifying framework that works for different load and generator models. Edges 
correspond to transmission lines. Each transmission line is characterized by its admittance Gij — jBij (j = y/—l). However, 
since we assume that transmission lines are lossless, we set Gij = 0 and thus each line is characterized by a single real and 
positive parameter B^j (the negated line susceptance). 
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Let Vi,9i,Pi and Qi denote voltage (magnitude), phase, active and reactive injection at the bus i respectively. Let pi = 
log (Vi). Let 9ij = 9i — Oj,pij = Pi — Pj. Buses are of three types: 

• (P, V) buses where active power injection and voltage are fixed, while voltage phase and reactive power adjust as 
conditions (e.g. power flows) are variables. The set of {P, V) buses is denoted by pv. 

• (P, Q) buses where active and reactive power injections are fixed, while voltage phase and magnitude are variables. The 
set of {P, Q) buses is denoted by pq. 

• Slack bus , a reference bus at which the voltage magnitude and phase are fixed, and the active and reactive power injections 
are free variables. The slack bus is denoted by S. 

0 and V (with no subscript) denote the vectors of phases at all buses except the slack bus and voltage magnitudes at all (P,Q) 
buses, p, V are indexed hy {i : i G pq}. For every bus i define Bi = ^ij- 

Let S be an ordered set of indices. We use the notation 



where i,j G S, to denote an [S'! x [S'! matrix whose rows and columns are indexed by the entries of S, with Mu = a,Mij = 
b, Mji = c, Mjj = d and all other entries equal to 0. Similarly, we use the notation M = [a]^ to denote an [S'! x [S'! matrix 
with the {i, i) — th entry equal to a and all other entries equal to 0. 

Given a vector x with indices i G S, denote by diag (x) the |S'| x [S'! matrix with diagonal entries Xi. Given a set C C M", 
Int (C) denotes the interior of the set. 

B. Background 

In this section, we introduce two power systems concepts that are critical for what follows: PF equations and energy functions. 

1) Power Flow Equations: PF Equations model the flow of power over the power system network. They are a set of coupled 
nonlinear equations that follow from Kirchoff’s laws applied to the AC power network. Circuit elements in the standard power 
systems models are all linear, if one ignores discrete elements like phase shifters and tap-changing transformers. Even though 
Ohm’s law and Kirchhoff’s law are linear in voltages and currents, power is their product and hence quadratic. Since we are 
interested in power, e.g., to determine power generation, specify demand, minimize line loss, PE equations are nonlinear. 

PE equations assume that the network is balanced, that is, the net sum of powers consumptions, injections and power 
dissipated is zero. This relies on the assumption that at the time-scale where the PE are solved (every 5 minutes or so), the 
system is in a quasi-steady state, i.e. the dynamic disturbances have been resolved through actions of the automatic control 
schemes represented by the automatic voltage regulators, power system stabilizers and primary and secondary frequency control 
systems. 

Since our main concern here is to study properties of the PE equations, we simply state them without derivation. At each 
node in the power system, there are four variables: voltage magnitude (Vi), voltage phase (6i), active power injection (Pi), 
reactive power injection (Qi). When solving the PE equations, at each node, two of these variables are fixed and the other two 
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are free. This leads to 2n equations in 2n variables, where n is the number of buses. Typically, power systems have two types 
of nodes: {P, V) nodes and {P, Q) nodes. At the {P, V) buses, the active power injection Pi and the voltage magnitude 
Vi = exp(/5i) are fixed. The fixed voltage magnitude is denoted Vi, to distinguish it from variable voltages at other nodes. 
Generators are typically modeled as (P, V) buses, since they have voltage regulators that inject reactive power to maintain 
voltages at fixed values. We can choose Vi = 1 for each i G pv, since any non-unit voltage set-point can be absorbed into the 
susceptances Bij in the active and reactive power flow equations. 

At the (P, Q) buses, the active power injection Pi and the reactive power injection Qi are held fixed. Loads are typically 
(P, Q) nodes, where Pi and Qi represent the active and reactive power demands (assumed constant). 

Usually, the variables Qi at the (P, V) buses are not considered explicitly. These variables can be obtained from the reactive 
power balance equations by plugging in there the values of {p,9). Finally, we have a special bus (slack bus) which absorbs 
all the power imbalances in the network. The injections Ps,Qs this bus are free variables whose values are adjusted to 
guarantee the power balance and 9s = Ps = 0. 

The power balance equations at all the buses, a set of coupled nonlinear equations in {p, 9), together constitute the PF 
equations: 


Pi = ^ Bij exp (pj -I- pj) sin (%), i e pv U pq 

jr^i 

(la) 

Pi = 0, i G pv 

(lb) 

Qi = ^ Bij (exp (2pi) - exp (pi -1- pj) cos (%)), i G pq 

(Ic) 

9s = 0, P5 = 0. 

(Id) 


2) Energy Functions: Energy functions were first introduced in the context of first integral analysis of the power system 
swing dynamic equations IIBH81IIICRP85I . The energy function consists of two terms - kinetic and potential. It is also a 
Lyapunov function for the swing dynamics. In this paper, we are not explicitly concerned with the dynamics, so the most 
useful property of the energy functions is that the stationary points of the energy function potential part map to solutions of 
the PF equations. The PF equations Q can be re-written in the following variational form 

dE 

([T^ : 0 = — Vi e pq U pv (2a) 

Uui 

dE 

( fT^ : 0 = -— Vi S pq (2b) 


Pi = 0 Vi e pv, P5 = 6*5 = 0 


(2c) 
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where E is called the energy function: 


E{p,9) = - ^ PiOi - ^ Q^pi 

zGpvUpq ^^Pq 

/exp {2pi) + exp {2pj) 


+ X! '^*•5 




“ X! 5iiexp(pi+pj)cos(%) 
{i,j)es 


Note that pi is a variable for i S pq and pi = 0 for z G pv, ps = 0, 6*5 = 0. 


(3) 


III. Energy Function and Convexity 

In this section, we study the convexity of the energy function ([^. We first discuss why convexity is an interesting and useful 
property, and then describe our results characterizing the domain of convexity of the energy function ([^. Before we proceed, 
let us define the domain of convexity precisely. 

Definition 1. The domain of convexity of the energy function is defined as: 


V = {{p,e) ■.V^E{ap,ae)>Q VaG[0,l]} 

Remark 1. In general, there could be several (disconnected) regions over which the energy function is convex IIASV81II . We 
require that for any (p, 9) G V, the energy function is convex at every point on the line segment connecting the origin to (p, 9). 
This ensures that we pick the region of convexity that includes the origin. 


A. Why is Convexity Interesting? 

The major contribution of this paper is the characterization of the convexity domain of the energy function (|^. We now 
discuss why this characterization is important and interesting. The principal reason is that solutions of the PF equations are 
stationary points of the energy function (|^. We will show, in section III-C that one can construct a convex optimization 
problem that will either find a solution of the PF equations within the convexity domain, or otherwise certify that there exist 
no solutions within the convexity domain. In this Section, we justify restricting ourselves to PF solutions contained in V, as 
the solutions in T) show some additional desirable properties. 

1) Asymptotic Stability: The dynamics of the power system, under certain reasonable assumptions, reduce to the so-called 
“swing dynamics”, see e.g. IIVWC85II . In this model, each generator is modeled as a constant voltage source behind a transient 
reactance. For each generator, there is an “internal” node (modeled as a constant voltage source subjected to the second order 
angular dynamics) and a “terminal” node (a node with 0 active and reactive power injections). We denote the set of internal 
generator nodes by pvg. The swing equation at an internal node i G pvg is 

dE{p,9) 


Mi9i + Di9i — — 


d9i 


(4) 


At all other nodes i G pv U pq, we have the algebraic active and reactive power balance equations Q. The dynamic state 
variables are {0^ : i G pvg} U {0^ : * G pvg}. 
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Theorem III.l. Let (p, 6) be an equilibrium of the swing dynamics such that (p, 9) G T> and satisfy all the algebraic equations 
Q- Then, (p, 9) is asymptotically stable, that is, linearizing 0 around (p, 9) results in a stable linear system. 

Proof. This follows from Lemma 5.1 in IITAV85II . which is a stronger result. □ 

2) Existence of Solutions: The following result is contained in our recent manuscript IDMI : 

Theorem III.2. Suppose that the network is a tree. Then, the PF equations 0 have a solution if and only if a solution is 
contained within T). 

Thus, at least for trees, it is sufficient to look for solutions of the PF equations contained in T). Preliminary numerical 
evidence (section [0 suggests that this observation also extends to mesh networks. Further, as we shall illustrate numerically in 
Section the solutions within T) are the “desirable” solutions (they have small phase differences between neighboring buses 
and sufficiently large voltage magnitudes at all (P, Q) buses). 


B. Characterization of Domain of Convexity of the Energy Function 

We now derive our main results on convexity of the energy functions. Since stationary points of the energy functions 
correspond to solutions of the PF equations, which are known to have multiple isolated solutions, the energy function is not 
globally convex. Thus, we focus on searching for restricted domains within which the energy functions are convex. In particular, 
we will characterize a subset of D such that phase differences between neighbors are smaller than 90 deg, since this is typical 
for normal power systems operations. It was shown in 0BBC13I that this restriction is actually sufficient for convexity in the 
case when all buses are (P, V) . However, when the network also contains (P, Q) buses additional conditions are required. 

Theorem III.3. The energy function E(p,9) is jointly convex in {p,9) over the convex domain CCD given by 


%- 

E 

iGpq 


9,1 <2 V(*,j)e£- 


2B, - yG 


,GpvU{5} 


exp {pji) 
' cos {9ji) 


- E 

(*j)6£,*,i6pq 


COS {9ij) 



exp(pij) 


P 0 


(5) 


( 6 ) 


Further, suppose that the network has a tree topology. Then, the domain of function convexity, D, equals C. Also, the energy 
function is strictly convex on Int (C). 


Proof. See appendix section VIII 


□ 


A number of Remarks are in order. 


Remark 2. Simple special cases 

When no (P, Q) nodes are connected to each other, the second term of the semi-definite inequality 0 is an empty sum. Thus, 















the semi-definite inequality is equivalent to 


2B, = 2[Y,B,A> Y. 


iGpvU{5},i~z 


Vjexp{-p,) 
cos {6i — 9j ) 


Since Bi = suffices that Viexp {—pj) < 2 cos {9i — Oj) for this to be true, which is satisfied if Vj > .7vi, 

~(^j \ < f ■ This is a fairly reasonable restriction in practice. A similar condition was discussed in IITAV85I . However, it was 
imposed only at the equilibrium point where the energy function reaches its minimum, as a sufficient condition to guarantee 
(small deviation) asymptotic stability. The asymptotic stability follows from our results as well, which guarantees that the 
equilibrium point lies within the domain of convexity of the energy function and it is hence asymptotically stable. 

Also, if there are no (P, Q) nodes, the matrix inequality condition in (|^ is empty and hence the convexity condition reduces 
to requiring that all phases differences (over existing lines) are smaller than ^ (see IIBBC13II 1. 

When (P, Q) -{P, V) connections are present, the semidefinite inequality does not simplify. However, in section we 
show numerically that the conditions imposed by the convexity domain are not very restrictive and most practical power flow 
solutions lie within the convexity domain. 


Remark 3. Conservatism 

For the tree networks, we have shown that the convexity conditions are necessary and sufficient (assuming that we are interested 
only in the domain with phase differences over all the lines smaller than |). For meshed networks, these conditions are sufficient 
but, possibly, not necessary. This is due to the fact that our convexity analysis has treated the edge variables Oij as independent 
for each edge (i, j) € £ while in reality the variables are coupled since Oij = 9i — 9j. Analysis of the gap is left for future 
work. We present some numerical results indicating that the gap is small in section |V] 


C. Solving Power Flow Equations via Convex Optimization 

Although energy function was developed initially as a tool to assess dynamic stability, the convexity results also yields 
consequences useful in the context of a number of other important power system applications. We choose to emphasize here 
one of these applications: Solution of the PF equations, leaving the discussion of other potential applications for future work. 

As we have already seen, the energy function technique provides a variational characterization of the PF equations (|^. 
Therefore, finding solution of the PF equations is equivalent to finding a stationary point of the energy function. 

Corollary 1. Let S C Int (C). Then, the power flow equations 0 have a solution (p*,9*) G S if and only if the following 
optimization problem has its (unique) optimal solution in S: 


min E Ip, 9) 


(7) 


Proof If {p*,9*) G S solves (0, it satisfies VpE {p*,9*) = gE [p*,9*) = 0. Thus, it also satisfies the KKT conditions 
for the optimization problem Q, forcing the respective Lagrange multipliers to be 0. Hence, the solution is optimal for Q. 

Conversely, suppose that (p*, 9*) G S C int (C) solves 0. Then, by complementary slackness, the Lagrange multipliers are 
0 and hence the KKT conditions reduce to VpE {p, 9) — 0, VgE (p, 9) = 0. Thus, (p*, 9*) G C also solves 0. 
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Further, by strict convexity of the energy function E {p, 6) on Int (C), there can be at most one stationary point of E in 
Int (C), and hence, the power flow solution, if it exists, is unique. □ 

IV. Related Work 

In this section, we discuss and contrast our work with related previous work. 


A. Principal Singular Surfaces 

The work closest in spirit to our approach is presented in the series of papers MTS72cl . IITS72bl . IITS72al . The authors of 
0TS72cL IITS72bl . IITS72all consider lossless networks where all buses are (P, V) buses. They define Singular surfaces, which 
are closed surfaces where det (V^P) = 0. The Principal Singular Surface (PSS) is the unique surface that encloses the origin 
6 = 0 and the Principal Region is the region enclosed by it. The principal region coincides with our definition of V when 
all nodes are {P,V) nodes, apart from the additional requirement that the phase differences are smaller than Our results 
characterizing C can be seen as constructing convex inner approximations to the principal region. 

This analysis was extended in 0ASV811 . where a number of counterintuitive properties of the set of solutions of the power 
flow equations were shown , even in the simple case when all buses are (P, V) buses. Based on these results, it was conjectured 
in this paper that the principal region is nonconvex in general. Further, there are examples of networks for which the only 
solutions to the power flow equations are unstable solutions. However, it was conjectured that if there is a stable solution, there 
must be one in the principal region. 

Extending this analysis and verifying the conjectures for our general setting with (P, Q) nodes is a direction of future work. 


B. Convex Relaxations of OFF 

The authors of the series of recent papers 0LTZ13I . IIZT13I . IILTZ13I . 0Lowl4bl . IILowl4al have studied convex relaxations 
of Optimal Power Flow (OPF) over tree networks. Various conditions were uncovered under which the convex relaxations 
of the OPF problems are exact, in the sense that an optimal solution of the original OPF problem can be recovered from 
the solution of the convex relaxation. In this work, we primarily deal with the PF equations, not OPF. However, for lossless 
networks our approach applies to arbitrary topologies and can potentially be extended to the OPF formulation, following the 
logic first time sketched in IIBBC13I : 


max min 
P,QeT (p,9)es 


Pidi - y] QiPi 

zGpvUpq *^Pq 


+ X! 

{ij)e£ 


exp {2pf + exp {2pj) 


~ Bij exp {pi + pj) cos {dij) 


-M{P,Q) 


where £t (P, Q) is a convex cost on the injections and T, S are convex operational constraint sets on the injections and voltages, 
respectively and A > 0 is a positive scaling factor. This is a convex-concave saddle point problem, as long as 5 C C and can 
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be solved efficiently using ideas from duality IIBV09I . If the constraints in S are not binding at the optimal solution, the inner 
minimization effectively enforces the PF equations Q, since the minimization corresponds to setting the gradient of the energy 
function to 0. The outer maximization then optimizes a cumulative cost consisting of the negative injection cost with an extra 
term added - energy function value at the optimum over line hows (corresponding to solution of the PF). We can choose A 
large enough so that the energy function term can be neglected for the outer maximization, so that we are effectively solving 
the OPF problem: 


min i (P, Q) 

p,Q,P,e 

Subject to 0, (P, Q) G T, (p, 9) G S 

Further developments on this line of work, along with conditions when our minimax formulation is equivalent to the OPF 
problem above, will be pursued in subsequent work. 

C. Conditions on Existence of solutions to PF equations 

Several papers have studied conditions for existence of solutions to the Power Flow Equations BBZ14IIIILSP99IIIWK82i . In 
IIMLD121 . the authors propose a sufficient condition for the insolvability of power how equations based on a convex relaxation. 
However, our approach differs from these in the following important ways: 

• To solve the PF equations in C, we provide necessary and sufficient conditions, that is, our approach finds a solution in 
C if and only if there exists a solution in C. 

• Our approach is algorithmic, that is, we provide an algorithm (based on convex optimization) that is guaranteed to find 
the solution efficiently (in polynomial time). 

• If there are additional operational constraints (p, 9) G S, with S C C, e.g. correspondent to line protection BSHOll and/or 
line how limits, we can additionally answer the question of whether there exists a PF solution with (p, 9) G S. This is 
an important contribution, since most of the time system operators are interested in finding power how solutions that 
additionally satisfy operational constraints. 

In IIBZ14II . the authors also propose an algorithm based on a contraction mapping consideration. However, the algorithm only 
works in a small ball around the origin in the {P, Q) space. Our results are stated in terms of a nonlinear convex constraint in 
(p, 9) space. Precisely understanding the set of {P, Q) for which the solution (p, 9) G C, and conversely the set of (p, 9) for 
which (P, Q) lies in a certain ball, is still an open problem, even for the special case where all buses are (P, V) buses. This 
setting was studied and the results were exteded and connected with results on synchronization in coupled oscillators in a series 
of recent papers IIDB12llflDCB13IIIDB14l . The authors provide distinct sufficient and necessary conditions on the injections for 
the existence of power how solutions with phase differences satisfying certain bounds. The authors also note that obtaining 
non-conservative sufficient conditions is still an open problem, and show numerically that a necessary condition is “almost” 
sufficient for most power systems. Extending this analysis to the setting with (P, Q) buses and obtaining non-conservative 
sufficient conditions on the injections that guarantee existence of power how solutions is an interesting direction for future 


work. 
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V. Numerical Illustrations 

This Section presents numerical experiments illustrating theoretical results resented above. 

A. Existence of Solutions: 2-bus Network 

Here we analyze the toy case of a two bus network with one generator and one load. The generator also serves as the slack 
bus. Thus, the only variables are the voltage phase 9 and voltage magnitude exp (p) at the load bus, and the parameters are the 
active power demand P and the reactive power demand Q at the load bus. The results, plotted in the figure]^ show the power 
flow solutions shown on the heat (iso-line) map of the energy function. There are two solutions for small loads. The convexity 
domain is the region inside the dashed green arc. The solution plotted as a green dot is a dynamically stable solution in the 
interior of the convexity domain. The solution plotted as a red dot is a dynamically unstable solution outside the convexity 
domain. As the loads are increased, the two solutions move closer to each other. At a certain critical value, the two solutions 
hit the boundary of the convexity domain and then vanish, thus showing that, beyond this critical loading level, there are no 
longer any solutions to the PF equations. 

This simple experiment supports our conjecture that in general, if there is a solution to the PF equations, there must be 
one found within the convexity domain. As stated earlier, this result was proven for networks with tree topology ODMI . In the 
special case when all nodes are the {P, V) buses, this result was conjectured (without proof) in 0ASV811 . with the caveat 
that one is only looking for dynamically stable solutions. We plan to investigate this conjecture for the general case with 
{P,Q) buses in further work. 

B. Theoretical versus Numerical Domain of Convexity 

In this section, we compare numerically these regions of the energy function and the reduced energy function convexity. In 
order to be able to visualize these regions, we restrict ourselves to a 3 bus toy system with one PV bus and two (P, Q) buses. 
Bus 1 is the PV bus (also the slack bus); its voltage magnitude Vi is fixed to 1 pu and its phase is selected as a reference, 
9i = 0. Buses 2 and 3 are {P,Q) buses with variable voltage magnitude V 2 = exp (^ 2 ), fa = exp(p 3 ) and voltage phase 
^ 2 ,^ 3 - ^^ 2)^3 are determined by solving the reactive PF equations 

Q2 = Bi2 {V2 — V 2 cos ( 6 * 2 )) + B23 [V2 — V2V3 cos (02 — 03 )) 

Qs = Si3 (Vi - V 3 cos ( 03 )) + B23 (Vi - V 2 V 3 cos (02 - 03 )) 

Active power injections are chosen as Pi = 3.7, P 2 = —1.7, P 3 = —2 (these numbers are in the normalized units) and the 
reactive power demands are <52 = —1-05,(53 = —1-24. The susceptances on the lines are B 12 = 26.88, P 13 = 26.88. We 
consider two cases for B 23 that lead to different network topologies: In the first case buses 2, 3 are not connected so B 23 = 0 
and in the second case we choose B 23 = 16.67. 

Given 0, we can solve the reactive power flow equations ( [T^ for p: We denote the solution by p (0). We study the convexity 

domain of the reduced energy function E {p (0), 0) so that we can plot the results (since the reduced energy function is only 

a function of two variables, ( 0 i, 02 ))- 
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Energy Function for p=q=.01 Energy Function for paqs.l 



(a) P = Q = .01 


(b) p = g = .1 



Energy Function for p=q=.21 



(c) p = g = .2 (d) p = g = .25 

Fig. 1; Solutions of the PF equations (red and green dotes) shown within the heat (iso-line) map of the energy function. The 
green dot is a dynamically stable solution inside the convexity domain, and the red dot is a dynamically unstable solution 
outside the convexity domain. Two distinct solutions are sufficiently far from each other at the initial (moderate) load. The 
two solutions get closer to each as the load increases. The solutions approach the boundary of the region of convexity from 
opposite sides and disappear once they hit the boundary. 


For different values of 02 , ^*3 G [~f i f ] > compute the value of the energy function E by solving the reactive PF equations 
for p 2 , Ps (using the Newton-Raphson iterative algorithm) and then plug the results into E {p, 9). For any value of 02 , 6 ^ 3 , we 
also estimate the Hessian VgP (p(0),0) by means of numerical differentiation. This gives us the exact region over which the 
reduced energy function, E {p{9), 0), is a convex function of 0. We compare this to the domain C of the energy function (and 
thus reduced energy function) convexity predicted by our theorems. Since the theorems only predict joint convexity of E {p, 0), 
we deduce reduced convexity by finding the set of 0 such that {p (0) ,6) G C. 
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Convexity Region; No (P,Q) buses connected 


Convexity Region: Normal Loading 


Convexity Region; High Loading 


Convexity Region: Over critical load 



(a) No (P,Q) connected 


(b) Regular Load 


(c) High Load 


(d) Loading over critical 


Fig. 2; Actual vs Estimated Domain of Convexity. Dark blue marks the region of non-convexity, brown represents the domain 
of predicted convexity (which coincides with actual convexity for all cases). One can see that as loading increases, the domain 
of convexity shrinks and then vanishes, indicating that there is no stable PF solution within the specified domain. 

Fig. 1^ shows the actual and predicted regions of the reduced energy function convexity. For the 3-bus case considered here, 
our theoretical results predict the region of convexity exactly. In Fig. 1^ we plot the convexity domain for the case B 23 = 0. 


All other Figures show the case of B 23 > 0 with a scaled-up injection. The scaling is 3 and 5.5 and 6 in Fig. 2b Fig. 2c and 


Fig. 1^ respectively. As the loading increases, the region of convexity shrinks and eventually we reach the critical load where 
the region of convexity becomes empty. At this point, there is no longer a stable power flow solution within the given domain, 
since the energy function must be convex around a stable equilibrium point. These results confirm that, at least for the 3-bus 
case, our theorems provide an exact description for the region of convexity. 


C. Power Flow In Larger Networks: Existence of Solutions 


In Section III-A2 we argued that one justification for seeking solutions only within the energy function’s domain of convexity 
is that there is an evidence to suggest that if there are no solutions within the Convexity Domain, then there are no solutions 
anywhere. We have a proof of this for the tree case IDMl and we conjecture that the result holds in general. In order to 
examine this hypothesis, we have studied PF solutions for two IEEE test networks; The IEEE 14 bus network and the IEEE 
118 bus network. We first modify these networks to become lossless, (we simply set conductances to 0 for all the transmission 
lines.) We start with the nominal base load profiles for these networks and gradually scale up the active power injections at 
all the buses (except the slack bus) by a factor k and scale up the reactive power injections by 6k at all the (P, Q) buses. 

In order to verify insolvability of the PF equations, we use the technique developed in 0MLD12I along with the code 
implementing this technique in the MATPOWER package OZMSTl II . For a fixed value of 6, we increase k until the technique 
from flMLD121 detects insolvability of the power flow equation. Our convex solver does not And a PF solution when the 
optimal solution to Q lies on the boundary of C or when the norm of the gradient of energy function with respect to (p, 6 ) is 
non-zero at the optimum. As k increases, we plot the norm of the gradient of the energy function at the optimum (a non-zero 
value indicates that there are no solutions in the convexity domain) and the value of the insolvability test (marking the result 
by 1 in the case of insolvability and 0 otherwise). 

For the IEEE 14 test case, the results are shown in Figures 3a|3b and 3c One observes that for 5 = 1, the scaling at which 
our convex approach fails to And a solution is exactly the point at which the technique from IMLD12I detects insolvability. 
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For smaller values of 6 , there is a slight gap between the two, but the values are still fairly close. 

For the IEEE 118 bus test case, the results are shown in Eigures 3a|3b and 3c We observe here (again) that the point at 
which 0MLD121 detects insolvability is close to the point at which our convex approach fails to find a power flow solution. The 
gap here is larger than the 14 bus case, but it is still rather small (difference of .1 in k, which amounts to a 2.5% difference 
in actual injections). 

The results overall show that our convex approach finds solutions to the power flow equations in almost all the cases when 
a solution exists. There is a potential gap close to the boundary of insolvability of the PE equations. We further observed in 
all the cases when our algorithm reported no solution, the Newton-Raphson based approaches implemented in MATPOWER 
also failed to And a solution. Based on this observation, we conjecture that the gap is actually due to the approximate nature 
of the 0MLD12I test (based on a sufficient but not necessary guarantees for the infeasibility). 


D. Operational Constraints and Convexity 

The convexity condition is a complicated nonlinear matrix inequality. In this Section, we describe the simple sufficient 

conditions for this inequality to be satisfied by computing bounds on the phase differences \9i — 9j and voltage magnitude 
ratios exp {\pi — Pj\) = max We Ax a bound on the log-voltage differences exp {\pi — Pj\) < bp and compute a 

bound on the phase differences \9i — 9 j \ < bg so that 

{{p,9) : exp(|pj - Pj\) < bp, \9i - 9j \ < bg V(i,j) S £} C C 


It is easy to see that the phase differences may be different for different lines, so we compute the tightest bound among all 
the lines. 


Eor the 14 bus case, the results are plotted in Eig. 5a The results show that if we allow voltage magnitude ratios across 
lines to be expdp^ — pj\) < 1.5, the phase differences can still be as large as 50deg. This is very reasonable for practical 
power systems, where these bounds are seldom exceeded. 


Eor the 118 bus case, the results are plotted in Eig. 5b Again, if we allow for voltage differences of up to exp {\pi — Pj\) < 
1.5, phase differences can still be as large as 45 deg. 

To conclude, our experimental results show that with fairly lax constraints on the voltage magnitude ratios and phase 
differences across neighboring buses, we are guaranteed to fall within the domain of convexity. In other words, most practical 
solutions to the PE equations lie within the domain of convexity. 


VI. Conclusions and Path Eorward 

This manuscript presents novel descriptions for the domain of convexity of the energy function of a lossless power 
transmission network. Since the energy function provides a variational description of the power flow equations, this allows us 
to develop an algorithm, based on a convex optimization, solving the PE equations. This enables various applications, including 
detecting if the PE equations have a solution in a certain domain. There is also a growing family of important power system 
applications which should benefit from these convexity and PE solvability results. These applications include (a) Optimal Power 
Plow, (b) State and Topology Estimation, (c) Distance to Insolvability/Pailure, (d) Optimal Load Shedding etc. 
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Another direction for future work is to extend these results to lossy networks. A special case, assuming that all lines within 
the network are of the same grade (the value of XijjRij = k W {i,j) G £), is discussed in Appendix 


VIII-B 


In the future 


work we plan to exploit continuity arguments to establish existence of solutions for networks with ^ values roughly constant. 
It should also be possible to quantify the degree of deviations that can be tolerated. 
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VIII. Appendix 


A. Proof of Theorem III.3 


Proof We introduce an independent edge variable 9k for each edge k G £. Let Og = {9k : k G £}. Let ki,k 2 denote the 
“from” and “to” ends of the edge k. We then have = A9. If the energy function is convex jointly in (p, 9, 9s), then it must 
be jointly convex in {p,9). Writing the energy function in terms of {p,9,9s), we derive E{p,9,9s) = ~ EiGpqUpv~ 
ZiGpq QiPi + Efe \Bk (exp (2pfci) + exp (2pfe2)) - Efe Bk exp {pki + Pk 2 ) cos {9k). In this expression only the first term 
depends on 9. Further, the first two terms depend linearly on (p, 9) and are hence convex. Therefore, convexity of the energy 
function reduces to convexity of the second and third terms. We study the Hessian which can be broken into 4 sub-matrices: 


/ VlE{p,9s) 


Kfi,E{p,9s) 

Vl,E{p,9s) 


L 

M 

jvj 

r 
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where M € Klpq|x|pql^ ]V e Rlpql^l^l, i? g For the matrix to be positive semi-definite, we require that R ^ 

0, M — NR~^N'^ y 0. It is easy to see that R is diagonal, with R^k = Bk exp {pki + Pk 2 ) cos {Ok). Further, 


Nik = \Aik\Bkeyip{pki -I-p/c^) sin(0^) G pq, fc G 

Aiij — ^ ^ I \Bk exp {pk^ Pk2 ) cos {Ok) j t 7 ^ J 

k 

Mu = 2Bi exp {2pi) - ^ \A,k\Bk exp {pk^ + Pfc^) cos {Ok) ■ 

k 

Let O = M - NR-^N^ G Mlpq|x|pql Then, 


for i ^ j and 


Oij — 


'^\A,k\\Ajk\Bkexp {pk^ + Pk 2 ) (cos {Ok) + 

k ^ 

'y \ Aik I|f?/c exp {pki + Pk 2 ) QQg (Qy,) 


sin^ {0k) \ 
cos {Ok) ) 


Note also that 


Ou = 2Bi exp {2pi) 


-^|^jfc|Bfeexp(pfc^ d-pfcs) 
k 


(^cos {Ok) + 


sin^ {0k) \ 
cos {Ok) j 


= 2Bi exp {2pi) -'^\Aik\Bk exp {pk^ + pfej 
k 


1 

cos {Ok) 




Lu = - \A,k\\Aik\Bk -^ 

exp(pi)exp(pj) 4^ cos {Ok) 


A 


La — 


O 


= 2B,-Y: M.*I exp in, + n, - 2ft) B, 


Then L ^ 0 if and only if O ^ 0 since L = diag (exp {—p)) Odiag (exp {—p)). Finally, note that L ^ 0 if and only if the 
following matrix is positive semidefinite; 

n pq 

2Bi— ^ Bij exp {pj — Pi) 


jepv,{i,j)es 


cos (0y ) 


E 


Bi 


ii,j)e£,i,jepq 


cos {Oi 


exp{pj - pA 


1 


exp(pi - Pj) 


To see that this is a convex constraint, we only need to observe that 


cos {Oij) 


exp {pj - Pi) 


exp {pi -pj) 


^3 


is ^-convex in {pj — Pi,0ij) (Lemma [^. The terms in the diagonal matrix are all concave functions of {p,0) and hence the 
rerms are ^-concave. Finally, plugging in Oij = Oi — Oj, one observes that the convexity is preserved and hence the energy 
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function is convex over C. Further, at (p, ff) = 0, the condition reduces to: 


M y 0, Mu = 2Bi — ^ Bij = ^ Bij,Mij = —Bij 

3~i 

Since M is symmetric and diagonally dominant, it must be positive semi-definite. Hence, 0 G C. Further, since C is convex 
{ap, aO) G C Va S [0,1]. Therefore, C CV. 

For the converse, we look at the tree networks. We assume that the tree is connected, else the energy function can be 
decomposed into a sum of independent functions on each connected component and the following arguments apply. First, note 
that for a connected tree with n buses, we find n—1 edges. Further, note that 9s = 0. Let 6 denote the vector of phases at all 
buses except at the slack bus. Then, there is a one-to-one correspondence between 6 and 9s, i.e., there exists an in— 1) x {n— 1) 
matrix A (submatrix of A formed by deleting the column corresponding to the slack bus) such that 9s = A9. Further, note 
that A is invertible. Hence, 9 = j 9s and the Hessian of E wrt (p, 9) is equal to 

( KEiP,9) yl^^^^^E{p,9s)A-^ \ 

[a-^\/1,^E{p, 9s f A-^Vl^^^E {p, 9,9s) A-^j 

Using Schur-complements and the invertibility of A, it is easy to see that the positive semi-definiteness of this matrix is 
equivalent to that of the matrix in the first part of the proof. The positive semi-definiteness of the bottom right block requires 
that \9s\ < If this inequality is strict, the block is positive definite and then the analysis in the first part is necessary and 
sufficient. The cases where this block is singular can be handled by a continuity argument. Hence, in this case, CAT), and 
combining the result from the first part, we arrive at I? = C. □ 


lemma 1. The matrix-valued function 


f[x,y) = 


exp {x) 1 

cos(y) I , / N 

' 1 exp y—x) 


is y-convex, that is / (Axi -f (1 - \)x 2 ,Xyi + (1 - \)y 2 ) A >^f{xi,yi) + (1 - A)/(a: 2 ,t/ 2 ) 


Proof, f is ^-convex if and only if tr (/ {x, y) R) is convex for all i? ^ 0 IIBV09II . Let R = 
2x2 positive semi-definite matrix. Then, we have 


^ 0 be an arbitrary 


b C ; 


tr(/(a;,y) R) = 


i exp (x) -\- c exp {—x) -1-26 
cos {y) 


aexp (a;)-I-cexp (—x) — 2 -y/oc b- 


ac 


cos {y) cos {y) 

Since R > Q, | 6 | < s/ac and hence the second term is convex (it is the inverse of a positive concave function, cos). The 

f 1 ■ (\Jae-x.-p(x/2) — ,/ce.-xp{ — x/2)\f . , , , \ r- / /n\ r- / ;r^\ i • 

hrst term can be rewritten as ^-cosfy)-check that ly'oexp (x/2) — y'cexp (—x/2) | is twice 

differentiable at all values of x and its second derivative is equal to 1-^0 exp (x/2) — -y/cexp (—x/2) |. Hence, this function is 
convex. Therefore, the function |i/aexp (x/2) — i/cexp (—x/2) p/cos(t/) is convex by the composition rules (since x^ jt is 
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convex in x,t, increasing in x and decreasing in t). Thus, / is ^-convex as long as \y\ < f • □ 

B. Lossy Networks 

Here we show that all of our results extend to the special class of lossy networks, where all the transmission lines have 

X G' 

the same/fixed value of the reactance-to-resistance ratio: = k. Further, we assume that all buses are (P, Q) buses, 

except for the slack bus. 

In this case, let us restate PF equations through the following linear combination of the PF equations written in their original 
form 


Pj + nQi = ^ + l) exp {pi + Pj) Bij sin (%), i G pq 




nPi - Qi = + l) exp {pi + Pj) Bij cos [Oij) Vi G pq 

Os = PS = 0 


(8a) 

(8b) 

(8c) 


Then the PF equations (|^ can be rewritten in a variational form: 


B{p,0) = ^ - ((Pi + nQi) 9i + [nPi - Qi) pi) 

iGpq 

+ V („> + 1) 

(i,i)€E ^ ^ ' 

- ^ Pij (K^ + l)exp(pi + pj)cos(6»i-6»j) 
(iV)6£ 


^p^ 


We arrive at the following statement. 


Theorem VIII.l. The energy function for the lossy system with constant -5 ratio (|9a|) is convex over the domain 


(9a) 

(9b) 


%-ef<- v(*,j)G£- 


E 

*Gpq 


OD Y- D exp(pj-,) 


i: 

(z,j)G£,ljGpq 


^ij 
COS {9i 


/ \ n P9 

exp (pji) 1 

1 exp(py)^ 


P 0 


( 10 ) 


( 11 ) 


Proof Almost identical to the proof of theorem III. 3 


□ 
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s=^ 



S=.5 



S=A 



Fig. 3: IEEE 14 Bus Case: Existence of Solutions as Load Scales. The solid black line is the norm of gradient of the energy 
function at the optimal solution. The value of k at which this becomes non-zero (marked by a dashed black line) is the point 
at which the convex solver fails. The red dashed line is the insolvability test. The value of k at which this becomes non-zero 
is the point at which insolvability is detected. 
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(5=1 



(a) 5 = 1 


(5=.5 



(b) 5 = .5 


(5=.1 



(c) <5 = .1 


Fig. 4: IEEE 118 Bus Case: Existence of Solutions as Load Scales. 















Convexity Region for IEEE 14 Bus System 

70 

60 

SO 

40 

30 

20 x 

10 

0 -'-'- 

1 1.5 2 2S 

Bound on exp(|p|-/’jl) 

(a) IEEE 14 Bus System 


Convexity Region for IEEE 118 Bus System 

60 - > ' ' " 


so 


4-40 


0 30 
■o 

I 20 

CO 

10 

ol-^-,-.-.- 

1 12 1.4 1£ 1.8 2 

Bound on exp(|/Y/>j|) 

(b) IEEE 118 Bus System 


Fig. 5; Convexity Domain and Operational Constraints 








